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We propose a new local algorithm for the thermalization of n- vector 
spin models, which can also be used in the numerical simulation of 
SU(N) lattice gauge theories. The algorithm combines heat-bath (HB) 
and micro-canonical updates in a single step — as opposed to the hy- 
brid overrelaxation method, which alternates between the two kinds 
of update steps — while preserving ergodicity. We test our proposed 
algorithm in the case of the one-dimensional 4-vector spin model and 
compare its performance with the standard HB algorithm and with 
other HB-inspired algorithms. 

I, INTRODUCTION 

The lattice formulation of quantum field theories provides a first principles ap- 
proach to investigate non-perturbative aspects of high-energy physics. In this for- 
mulation the theory becomes equivalent to a statistical mechanical model, which 
can be studied numerically using Monte Carlo simulations (see for example [1,2] 
and references therein). As a consequence, the system considered evolves according 
to a Markov process in the so-called Monte Carlo time and the action-weighted 
configuration-space average of the observables is substituted by a time average over 
successive (independent) field configurations of the system. The possibility to use 
Monte Carlo simulations to study the theory nonperturbatively is especially impor- 
tant in the case of quantum chromodynamics, the nonabelian SU(3) gauge theory 
describing the strong interactions between hadrons. These simulations are compu- 
tationally very demanding and must be done using local thermalization algorithms, 
since global methods (such as cluster algorithms) do not work well in this case. 

We should notice that, based on the above-mentioned equivalence, when the 
continuum limit of the lattice quantum field theory is taken, the corresponding sta- 
tistical mechanical model approaches its critical point. More precisely, the distance 
from the critical point is given by the lattice parameter (3 (directly related to the 
lattice bare coupling constant), which corresponds to an inverse temperature. For 
nonabelian gauge theories (as a consequence of asymptotic freedom) the continuum 
limit is given by (3 — > oo, i.e. the critical point corresponds to temperature zero. 
The expected critical behavior is therefore similar to the one of a two-dimensional 
continuous-spin model (e.g. the classical Heisenberg or n-vector spin model with 
n > 2) or of a one-dimensional spin model, which show a critical point only at 
temperature zero, or (3 — > oo. 

The process of obtaining independent field configurations is called thermalization 
and is usually carried out by applying at each link of the lattice a local algorithm, 
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such as Metropolis or heat bath (HB). When a critical point is approached, this pro- 
cess is afflicted by the well-known phenomenon of critical slowing- down (CSD) [3], 
which increases the correlation among successive field configurations. This implies 
that the integrated auto- correlation time Tint increases as a power of the lattice side 
N. In particular, for the Metropolis or HB algorithms one has Tint ~ N 2 , i.e. the 
dynamic critical exponent z is equal to 2. Since statistical Monte-Carlo errors are 
proportional to v / 2r int , numerical simulations become increasingly inefficient close 
to a critical point. In order to reduce the problem of CSD one can combine the stan- 
dard Metroplis and HB algorithms with so-called micro-canonical updates, allowing 
larger jumps in the configuration space and therefore improving the generation of 
independent samples. This is the idea behind the so-called hybrid overrelaxation 
method [4]. In general, adding a few micro-canonical sweeps greatly reduces CSD 
and, correspondly, the computational work. 

A modification of the heat-bath algorithm, called overheat-bath (OHB), was in- 
troduced some years ago in Ref. [5]. The basic idea was to incorporate a micro- 
canonical move directly into the heat-bath step, thus reducing the computational 
cost while preserving the large moves in the configuration space. As it turns out, 
combining the two moves (heat-bath and micro-canonical) in a single step leads to 
a significant improvement in performance when compared to the hybrid overrelax- 
ation method described above, which is based on alternating the two kinds of up- 
date moves. The resulting algorithm was indeed able to speed up the thermalization 
process, but, as already stressed in Ref. [5], it is not clear if it preserves ergodicity 
(especially when working at small temperatures). The OHB algorithm is used to- 
day in numerical simulations [6], usually combined with other algorithms in order 
to ensure ergodicity. In this work we propose a modification of the overheat-bath 
algorithm, which we call the modified heat-bath algorithm (MHB), incorporating a 
micro-canonical move into the heat-bath step without compromising ergodicity. In 
order to test the MHB algorithm — and compare its performance with the stan- 
dard HB algorithm and the OHB algorithm — we consider the 4- vector spin model 
on a 1-dimensional lattice. Note that the model is exactly solvable, which makes 
it especially suited for testing the algorithm and comparing results to the exact 
solution. 

Let us also note that, due to the isometry between the groups 0(4) and SU(2), 
it is possible to study the SU (2) case with the same algorithm used for the 4- vector 
case. More precisely, the local update — i.e. the update of a single spin or gauge 
field variable while holding the rest of the lattice fixed — of the SU(2) lattice gauge 
theory is identical to the one for the 4-vector model. This does not mean, however, 
that the global update of an SU (2) gauge field configuration can be obtained from 
the corresponding update of a 4-vector model. Indeed, the latter update can be 
done very efficiently using the Swendsen- Wang- Wolff algorithm, while no such class 
of algorithms works well for lattice gauge theories. The method proposed here is 
therefore intended mostly for use in simulations of lattice gauge theories and not of 
n- vector models, although hybrid overrelaxation methods have also been applied in 
large-scale simulations of n- vector models [7]. Of course, an efficient thermalization 
in the SU(2) case is of great physical interest, since the quenched SU(N C ) case (for 
N c > 3) is usually studied applying the S'[/(2)-embedding technique introduced in 
Ref. [8]. The SU(2) embedding is also used for simulating other Lie groups, such 
as the Sp(2) and Sp(3) groups, in studies of the deconfinement phase transition [9]. 
Preliminary results for the 2d SU(2) case have been presented in [10]. 
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II. THE 4- VECTOR SPIN MODEL AND THE ALGORITHMS 



The 4-vector spin model (on a 1-d lattice) is defined by the Hamiltonian 

N 

n = -(3 ■ s x+ i , (i) 

x=l 

where the spins S x are four-dimensional unit vectors, (3 ~1/Temperature , N is 
the lattice side and • indicates a scalar product. 

In the case of a local algorithm, one has to consider the contribution to the 
Hamiltonian due to a single spin S x . This gives H ss = —[3S X ■ H x + constant, 
where the "effective magnetic field" H x is given by H x = S x -i + S x +i . (Here we 
consider periodic boundary conditions.) The above equation can also be written as 

Hss = ~^~Y- Tr S x Hi + constant , (2) 

where S x and H x are now interpreted as SU(2) matrices in the fundamental 
representation and N x — V det H x . Note that Eq. (2) is clearly analogous to the 
expression of the single-link action obtained by considering the contribution of a 
single link variable to the SU(2) Wilson action [10]. 

Using the single-side action (2) and the invariance of the group measure under 
group multiplication, one obtains the HB update [11] 

S° ld -» S^ ew = UH X . (3) 

Here, U = uq! + i J2j u j <T j ls an SU{2) matrix, / is the 2x2 identity matrix, 
<jj are the three Pauli matrices and uo is randomly generated according to the 
distribution 

\Jl-ul exp ( f3 N x u ) du . (4) 

At the same time, the vector u — with components Uj normalized to yl — Uq — 
points along a uniformly chosen random direction in three-dimensional space [2]. 

In the hybrid over-relaxed algorithm one does m micro-canonical (or energy- 
conserving) update sweeps over the lattice, followed by one HB sweep. The micro- 
canonical update is a local deterministic transformation applied to the SU (2) matrix 
S x , which does not change the value of the Hamiltonian. This is achieved by 
considering the update 

- S°J d . (5) 

As stressed in the Introduction, this update represents a large move in the configu- 
ration space, allowing the hybrid over-relaxed algorithm to reduce CSD at the price 
of a greater computational cost, due to the micro-canonical sweeps. Usually, by 
setting m = 2, 3 one obtains a strong reduction in the value of Ti nt while increasing 
the computational cost by a factor smaller than 2. 

In the OHB [5] one tries to include the micro-canonical step (5) directly into the 
heat-bath algorithm. To this end one generates uo according to the distribution (4) 
while the components of u are not randomly chosen but are set using the relation 
u = —w, where W = wol + i<r ■ w = S° ld H x . As a final step, the vector u 
is normalized to \J\ — Uq. Clearly, the idea here is the same one applied in the 
standard hybrid over-relaxed algorithm: one tries to maximize the move in the 
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configuration space by changing the sign of the component of W that is orthogonal 
to the effective magnetic field. (Note that the action S — — /3 N/2TrW can be 
viewed as the action of a matrix W in an effective magnetic field given by the identity 
matrix I.) Clearly, this step does not obey the uniform distribution for u but is a 
microcanonical move. Indeed we can think of the algorithm as a two-step process: a 
HB move followed by a microcanonical step. Thus, the algorithm is exact but may 
not be ergodic. We analyze the conditions for applicability of the OHB algorithm in 
a separate work [12], but it is clear that for some initial configurations the algorithm 
can get "trapped" inside a subset of the space of configurations, compromising the 
ergodicity condition. One such configuration is, for example, a "cold" start for an n- 
vector model, in which all spins are aligned along a fixed direction. In this case, the 
OHB is not able to change the lattice configuration at all. There is a clear problem 
also if the initial configuration in the SU (2)-lattice-gauge-theory case is given by 
variables in an Abelian subgroup. As can be easily seen, the update will change 
the initial configuration in this case, but the resulting Markov chain will remain 
restricted to the Abelian sector, without exploring the full space of configurations. 

In this work we propose a modified HB algorithm (MHB) in which the generation 
of the SU{2) matrix U is done as in the HB case, but followed by one final step: 
if the scalar product u ■ w is positive then one does u — ► —u, where the matrix 
W has been defined above. This may also be thought of as a modification of the 
overheat-bath (OHB). As said above, the basic idea, in both cases, is to incorporate 
a micro-canonical move into the heat-bath step. The difference is that in our case 
the direction of u is randomly chosen and only its sign is modified, according to the 
above rule, while in the OHB algorithm one fixes u cx — w directly. Our modification 
should ensure ergodicity while keeping the advantage of having a micro-canonical 
step included into the HB step. We have indeed verified that the MHB algorithm 
has no problem with cold starts and shows a better r,„t than that of the OHB 
algorithm for energy-related observables. On the other hand, in the OHB case, 
the move in configuration space is more optimized and the iteration cost is lower, 
since the algorithm is simpler and it needs fewer random numbers. Our results are 
presented in the next section. 

III. RESULTS AND CONCLUSIONS 

In order to study the CSD of the new algorithm we have to investigate if, and with 
what exponent, its relaxation time Tint (for certain quantities we are interested in) 
diverges as the lattice size N increases. To this end, we have to evaluate Tint for 
different lattice sides N at "constant physics", namely as the lattice size is increased, 
the physical size of the lattice should remain fixed. This is done by introducing a 
correlation length £ and by keeping the ratio £/N fixed, with N <C £ in order 
to keep finite-size effects under control. For the 1-d 4-vector spin model one has 
£ cx (3 [13], thus one should tune N cx (3 in order to keep the ratio £/N constant. 
In our simulations we considered the lattice sizes N — 32, 64, 96, 128, 160, 224, 256 
and (3 = 2.5 N/32 and (3 = 5.0 TV/32. In all cases we did 1.1 x 10 6 thermalization 
sweeps and for the analysis we discarded the first 10 5 sweeps. 

In order to test the new algorithm we have evaluated the energy density e = 
N^ 1 J2x=i ' S x+ i , the magnetization M 2 = (J2x=i S x ) 2 an d the magnetic 
susceptibility \ — A -1 (At 2 ) . Let us recall that there exists an exact solution [13] 
for the 1-d 4-vector spin model, even at finite lattice side N and with periodic 
boundary conditions. Thus, we can check the numerical results for the various 
quantities against the exact solution (see also [14, Table 1]). 

In all cases we have evaluated the integrated auto-correlation time Tint using an 
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automatic windowing procedure [3] with two different window factors (c = 4 and 
8). We also employ a method based on a comparison between the naive statistical 
error with a binning error [15]. We checked that these three estimates are always 
in agreement. The error on the integrated auto-correlation time n n t has been 
evaluated using the Madras-Sokal formula [3]. 

Results are reported in Figs. 1 and 2. The data obtained for Tint have been fitted 
to the Ansatz T int — aN z in order to estimate the dynamic critical exponent z. 
Our study shows that the MHB algorithm has essentially the same critical exponent 
z of the standard HB algorithm but the value of the integrated auto-correlation 
time Tint is about 30% smaller compared to the standard HB algorithm. This 
implies a reduction in the statistical error and in the computational cost by a factor 
of about 20%. Similar results have been obtained in the SU{2) case [10]. Thus, 
a good decorrelation can be reached without the use of multiple micro-canonical 
sweeps. In particular, from our simulations it appears that the HB algorithm with 
to micro-canonical steps is essentially equivalent to the MHB algorithm with m — 1 
micro-canonical steps. This feature may be useful in parallel computing [2] due to 
the reduction of the inter-node communication. 

At the same time, from our data one sees that the MHB algorithm has a critical 
exponent z larger than that of the OHB algorithm. In order to reduce CSD for 
the MHB algorithm, while keeping its ergodicity properties, we are now testing an 
algorithm that interpolates between MHB and OHB. 

A more extensive study of these algorithms, both in the 0(4) and in the SU(2) 
cases, will be presented in a forthcoming work [12]. 
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FIG. 1. Integrated auto-correlation time of the susceptibility as a function of the lattice 
side N for HB (x), OHB (o) and MHB (o), with (left) and 1 (right) overrelaxation step. 
Data have been fitted to the Ansatz Ti„t = aN z . 
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FIG. 2. Integrated auto-correlation time of the susceptibility as a function of the lattice 
side N for HB (x), OHB (o) and MHB (o), with 2 (left) and 3 (right) overrelaxation step. 
Data have been fitted to the Ansatz Ti„t = aN z . 
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